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Abstract 

The purpose of this text is to give a self-contained description of the basic theory of the 
projector augmented-wave (PAW) method, as well as most of the details required to make 
the method work in practice. These two topics are covered in the first two sections, while 
the last is dedicated to examples of how to apply the PAW transformation when extracting 
non-standard quantities from a density-functional theory (DFT) calculation. 

The formalism is based on BlochPs original formulation of PAW [1] , and the notation and 
extensions follow those used and implemented in the GPAw[2] code. 
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1 FORMALISM 



1 Formalism 

By the requirement of orthogonality, DFT wave functions have very sharp features close to the 
nuclei, as all the states are non-zero in this region. Further out only the valence states are non-zero, 
resulting in much smoother wave functions in this interstitial region. The oscillatory behavior 
in the core regions, requires a very large set of plane waves, or equivalently a very fine grid, to 
be described correctly. One way of solving this problem is the use of pseudopotentials in which 
the collective system of nuclei and core electrons are described by an effective, much smoother, 
potential. The KS equations are then solved for the valence electrons only. The pseudopotentials 
are constructed such that the correct scattering potential is obtained beyond a certain radius 
from the core. This method reduces the number of wave functions to be calculated, since the 
pseudo potentials only have to be calculated and tabulated once for each type of atom, so that 
only calculations on the valence states are needed. It justifies the neglect of relativistic effects in 
the KS equations, since the valence electrons are non-relativistic (the pseudopotentials describing 
core states are of course constructed with full consideration of relativistic effects). The technique 
also removes the unwanted singular behavior of the ionic potential at the lattice points. 

The drawback of the method is that all information on the full wave function close to the nuclei 
is lost. This can influence the calculation of certain properties, such as hyperfine parameters, and 
electric field gradients. Another problem is that one has no before hand knowledge of when the 
approximation yields reliable results. 

A different approach is the augmented- plane- wave method (APW), in which space is divided into 
atom-centered augmentation spheres inside which the wave functions are taken as some atom-like 
partial waves, and a bonding region outside the spheres, where some envelope functions are defined. 
The partial waves and envelope functions are then matched at the boundaries of the spheres. 

A more general approach is the projector augmented wave method (PAW) presented here, which 
offers APW as a special case[3], and the pseudopotential method as a well defined approximation [4]. 
The PAW method was first proposed by Blochl in 1994[1]. 

1.1 The Transformation Operator 

The features of the wave functions, are very different in different regions of space. In the bonding 
region it is smooth, but near the nuclei it displays rapid oscillations, which are very demanding 
on the numerical representation of the wave functions. To address this problem, we seek a linear 
transformation T which takes us from an auxiliary smooth wave function \ip n ) to the true all 
electron Kohn-Sham single particle wave function | ip n ) 

\i>n)=t\i>n), (1) 

where n is the quantum state label, containing a k index, a band index, and a spin index. 
This transformation yields the transformed KS equations 

-pHf\i>„) = ejifltn), (2) 

which needs to be solved instead of the usual KS equation. Now we need to define T in a suitable 
way, such that the auxiliary wave functions obtained from solving (2) becomes smooth. 

Since the true wave functions are already smooth at a certain minimum distance from the core, 
T should only modify the wave function close to the nuclei. We thus define 

t=l + £f°, (3) 

a 

where a is an atom index, and the atom-centered transformation, T a , has no effect outside a certain 
atom-specific augmentation region |r — R a | < r". The cut-off radii, r" should be chosen such that 
there is no overlap of the augmentation spheres. 



1.1 The Transformation Operator 



3 



Inside the augmentation spheres, we expand the true wave function in the partial waves (f>f, 
and for each of these partial waves, we define a corresponding auxiliary smooth partial wave 4>f, 
and require that 

ic> = (i+t a M) & t a \tf) = m-\fc) (4) 

for all i,a. This completely defines T, given <f> and 4>. 

Since T a should do nothing outside the augmentation sphere, we see from (4) that we must 
require the partial wave and its smooth counterpart to be identical outside the augmentation sphere 

Va, #(r) = #(r), for r > r a c 

where <pf(r) = (r|0?) and likewise for <j>f. 

If the smooth partial waves form a complete set inside the augmentation sphere, we can formally 
expand the smooth all electron wave functions as 

|^n) = E P «il^)' fOT l r - Ra |< r 'c ( 5 ) 

i 

where are some, for now, undetermined expansion coefficients. 
Since \<pf) — T|0") we see that the expansion 

\^ n ) = f\i> n ) = KM), for |r - R a \ < r a c (6) 

i 

has identical expansion coefficients, P^. 

As we require T to be linear, the coefficients P^ must be linear functionals of \ipn), i-e. 

Ki = m\$n) = J dvpT{v - R a )Mr), (7) 

where \pf) are some fixed functions termed smooth projector functions. 

As there is no overlap between the augmentation spheres, we expect the one center expansion 
of the smooth all electron wave function, |^°) = J2i \&i)(Pi IV'n) to reduce to \tp n ) itself inside the 
augmentation sphere defined by a. Thus, the smooth projector functions must satisfy 

i 

inside each augmentation sphere. This also implies that 

(ptM) = Si ui , ,for|r-R°|<r c ° (9) 

i.e. the projector functions should be orthonormal to the smooth partial waves inside the augmen- 
tation sphere. There are no restrictions on outside the augmentation spheres, so for convenience 
we might as well define them as local functions, i.e. pf(r) — for r > r°. 

Note that the completeness relation (8) is equivalent to the requirement that pf should produce 
the correct expansion coefficients of (5)-(6), while (9) is merely an implication of this restriction. 
Translating (8) to an explicit restriction on the projector functions is a rather involved procedure, 
but according to Blochl, [1], the most general form of the projector functions is: 

(p7i = E({(/»« 8 >})«(/fi. ( 10 ) 

3 

where is any set of linearly independent functions. The projector functions will be localized if 
the functions |/") are localized. 

Using the completeness relation (8), we see that 



f a = Y,T a \ti)(p? i = E G#> - i#»<p?i. 

i i 
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where the first equality is true in all of space, since (8) holds inside the augmentation spheres 
and outside T a is zero, so anything can be multiplied with it. The second equality is due to (4) 
(remember that |0°) — \<f>f) — outside the augmentation sphere). Thus we conclude that 

t=i+x;E(i#>-i#»<i5?i- (ii) 

a i 

To summarize, we obtain the all electron KS wave function ip n (r) = (r\ijjn) from the transfor- 
mation 

Mr) = Mr) +EE (#W - (I 2 ) 

a i 

where the smooth (and thereby numerically convenient) auxiliary wave function ipn(r) is obtained 
by solving the eigenvalue equation (2). 

The transformation (12) is expressed in terms of the three components: a) the partial waves 
<f>i(r), b) the smooth partial waves </>"(r), and c) the smooth projector functions pf(r). 

The restriction on the choice of these sets of functions are: a) Since the partial- and smooth 
partial wave functions are used to expand the all electron wave functions, i.e. are used as atom- 
specific basis sets, these must be complete (inside the augmentation spheres), b) the smooth 
projector functions must satisfy (8), i.e. be constructed according to (10). All remaining degrees of 
freedom are used to make the expansions converge as fast as possible, and to make the functions 
termed 'smooth', as smooth as possible. For a specific choice of these sets of functions, see section 
2.2. As the partial- and smooth partial waves are merely used as basis sets they can be chosen as 
real functions (any imaginary parts of the functions they expand, are then introduced through the 
complex expansion coefficients P^)- In the remainder of this document <fi and <fi wm t> e assumed 
real. 

Note that the sets of functions needed to define the transformation are system independent, 
and as such they can conveniently be pre-calculated and tabulated for each element of the periodic 
table. 

For future convenience, we also define the one center expansions 

^(r) = X>?(r)<p^„>, (13a) 

i 

C(r)=£^(r)(^ n >. (13b) 

i 

In terms of these, the all electron KS wave function is 

Mr) = Mr) + «( r - RQ ) - C( r - RQ ))' (I 4 ) 

a 

So what have we achieved by this transformation? The trouble of the original KS wave functions, 
was that they displayed rapid oscillations in some parts of space, and smooth behavior in other 
parts of space. By the decomposition (12) we have separated the original wave functions into 
auxiliary wave functions which are smooth everywhere and a contribution which contains rapid 
oscillations, but only contributes in certain, small, areas of space. This decomposition is shown on 
the front page for the hydrogen molecule. Having separated the different types of waves, these can 
be treated individually. The localized atom centered parts, are indicated by a superscript a, and 
can efficiently be represented on atom centered radial grids. Smooth functions are indicated by a 
tilde ~. The delocalized parts (no superscript a) are all smooth, and can thus be represented on 
coarse Fourier- or real space grids. 

1.2 The Frozen Core Approximation 

In the frozen core approximation, it is assumed that the core states are naturally localized within 
the augmentation spheres, and that the core states of the isolated atoms are not changed by the 
formation of molecules or solids. Thus the core KS states are identical to the atomic core states 
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K) = K' corc >, 

where the index n on the left hand site refers to both a 
specific atom, a, and an atomic state, a. 

In this approximation only valence states are included 
in the expansions of \ipn), (6), and \ip n )i (5). 

Figure 1, shows the atomic states of Platinum in its 
ground state, obtained with an atomic DFT program 
at an LDA level, using spherical averaging, i.e. a spin- 
compensated calculation, assuming the degenerate occu- 
pation 9/10 of all 5d states, and both of the 6s states 
half filled. It is seen that at the typical length of atomic 
interaction (the indicated cut-off r c = 2.5 Bohr is approx- 
imately half the inter-atomic distance in bulk Pt) , only 
the 5d and 6s states are non-zero. 

1.3 Expectation Values 

The expectation value of an operator O is, within the frozen core approximation, given by 

val core 

(d) = £/„<woiVn> +EE^' coro i6iC' corc >- ( 15 ) 

n a a 

Using that (ip n \0\tfj n ) = (ijj n \T' OT\ijj n ) , and skipping the state index for notational convenience, 
we see that 

(fl>\d\i>) = (i> + E(V>° - i> a )\6\i> + Ew»° - ^ a )> 

a a 

= (4>\d\^) + - r\o\r' - ¥) + E (&\o\r 

aa' a 

= $|3|$ +e ((r\d\r) - md\r)) 

a 

+ E (V - ^ a i^ - ^ + ^ Q io^ a 

a 

For local operators 1 the last two lines does not contribute. The first line, because \ip a — ip a ) is only 
non-zero inside the spheres, while \ip — ip a ) is only non-zero outside the spheres. The second line 
simply because \tp a — ip a ) is zero outside the spheres, so two such states centered on different nuclei 
have no overlap (provided that the augmentation spheres do not overlap). 
Reintroducing the partial waves in the one-center expansions, we see that 

val val val 

E fnmo\r„) = E /» E^a mca) = X>£ E upzik^ m 

and likewise for the smooth waves. 

Introducing the Hermitian one-center density matrix 

Ki, = E /« P "n P ™ 2 = E /« K ) (Pi l^n)- (18) 



Local operator O: An operator which does not correlate separate parts of space, i.e. (r|0|r'} = if r ^ r'. 



— Is 




"DO 0.5 1.0 1.5 2.0 2.5 3.0 3.5 
t [Bohr] 



Figure 1: The core states of Platinum 



$») + (ip a -ip a \0\ 



(16) 
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We conclude that for any local operator O, the expectation value is 

val 

fn$n\0\i> n ) +] 



(6) =£/n<4Wn>+ EE («i%£> - (<|6fe>) z?? li2 + EE(^ coro ioic core >. 

a a 

(19) 



1.4 Densities 

The electron density is obviously a very important quantity in DFT, as all observables in principle 
are calculated as functionals of the density. In reality the kinetic energy is calculated as a functional 
of the orbitals, and some specific exchange-correlation functionals also rely on KS-orbitals rather 
then the density for their evaluation, but these are still implicit functionals of the density. 

To obtain the electron density we need to determine the expectation value of the real-space 
projection operator |r)(r| 

n(r) =E/n<^n|r)<r|V>„) = E/«l^n(r)| 2 , (20) 

n n 

where /„ are the occupation numbers. 

As the real-space projection operator is obviously a local operator, we can use the results (19) 
of the previous section, and immediately arrive at 

val core 

n(r) - e /«i^i 2 + E E (« - cc) D u, + E E i^ ,core i 2 - ( 21 ) 

To ensure that (21) reproduce the correct density even though some of the core states are 
not strictly localized within the augmentation spheres, a smooth core density, h c (r), is usually 
constructed, which is identical to the core density outside the augmentation sphere, and a smooth 
continuation inside. Thus the density is typically evaluated as 

n(r) = n(r) + EK«-^ a (r)), (22) 

a 

where 

val 

n(r) = E/«l^( r )| 2 + ^( r ) ( 23a ) 

n 

n%r) = E^CWCM +<W ( 23b ) 

M2 

n a (r) =J2 D liA(*)ti( r ) +«cW (23c) 



1.5 Total Energies 

The total energy of the electronic system is given by: 

E[n] = T.[n] + U H [n] + V ext [n] + E xc [n}. (24) 

In this section, the usual energy expression above, is sought re-expressed in terms of the PAW 
quantities: the smooth waves and the auxiliary partial waves. 

For the local and semi-local functionals, we can utilize (19), while the nonlocal parts needs 
more careful consideration. 
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1.5.1 The Semi-local Contributions 

The kinetic energy functional T s = J2 n /n(VVi| — 5 V 2 |?/> n ) is obviously a (semi-) local functional, so 
we can apply (19) and immediately arrive at: 

T s [{il n }]=Y J Ui>n\-\V 2 \^n) 

n 

val (25) 

= E fn$n\ \^n) + E ( T c + AT ^ ll2 ) > 
n a 

where 

core 

Tc = E^" 00 "! - lV 2 |C' cor °> and AT? li2 = « | - ±V 2 |C 2 > - (< I - ^V 2 |C 2 ). (26) 

a 

For LDA and GGA type exchange-correlation functionals, E xc is likewise, per definition, a semi-local 
functional, such that it can be expressed as 

E xc [n] = E xc [n] + E (E xc [n a ] - E xc [n a ]) . (27) 
a 

By virtue of (23b)-(23c) we can consider the atomic corrections as functionals of the density matrix 
defined in (18), i.e. 

E xc [n] = E xc [n] + E A^ C [{^,J], (28) 

a 

where 

^E a xc [{Dl i2 }\ = E xc [n a ] ~ E xc [n a ]- (29) 



1.5.2 The Nonlocal Contributions 

The Hartree term is both nonlinear and nonlocal, so more care needs to be taken when introducing 
the PAW transformation for this expression. 

In the following we will assume that there is no 'true' external field, such that 14 x t[n] is only 
due to the static nuclei, i.e. it is a sum of the classical interaction of the electron density with the 
static ionic potential, and the electrostatic energy of the nuclei. 

We define the total classical electrostatic energy functional as 

E c [n] - U H [n] + V cxt [n] = I((n)) + (n| £ a Z a ) + \ £ Q#a ^ a |^ a '), (30) 
where the notation (f|g) indicates the Coulomb integral 

and I have introduced the short hand notation ((/)) = (f\f)- In (30), Z a (r) is the charge density 
of the nucleus at atomic site a, which in the classical point approximation is given by 

Z a (r) = -Z a S(r-R a ) (32) 

with Z a being the atomic number of the nuclei. As the Hartree energy of a density with non-zero 
total charge is numerically inconvenient, we introduce the charge neutral total density 



p(r) = n(r) + E ^ ( r ) (= "electrons + "-nuclei)- 
a 



(33) 
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In terms of this, the coulombic energy of the system can be expressed by 

E c [n]=U' H [p] = \{{n + Y, a Z a ))' (34) 

where the prime indicates that one should remember the self-interaction error of the nuclei 
introduced in the Hartree energy of the total density. This correction is obviously ill defined, and 
different schemes exist for making this correction. As it turns out, this correction is handled very 
naturally in the PAW formalism. 

For now, we will focus on the term ((/?)) = ((n + E a Z a )). If one where to directly include the 
expansion of n(r) according to (22), one would get: 

((n + E a Z a )) = ((n + Ea "° - « a + za )) 

= ((»)) + - " a + za \ n<1 ' - ™ a ' + za ') + 2 X^i™ Q - " a + z<1 )' 

aa' a 

where in the last expression, the first term is the Hartree energy of the smooth electron density, 
which is numerically problematic because of the nonzero total charge. The second term contains a 
double summation over all nuclei, which would scale badly with system size, and the last term 
involves integrations of densities represented on incompatible grids (remember that the one-center 
densities are represented on radial grids to capture the oscillatory behavior near the nuclei) 2 . This 
is clearly not a feasible procedure. To correct these problem we add and subtract some atom 
centered compensation charges Z a : 

((n+ E a Z a + E a [Z a - Z a ] )) = {{n+ E a Z a )) + Eaa' K - n a + Z a - Z a \n a ' - n a ' + Z a ' - Z a ) 

+ 2 + Ea' za '\ na -n a + Z a - Z a ). 

a 

If we define Z a (r) in such a way that n a (r) — h a (r) + Z a (r) — Z a (r) has no multipole moments, i.e. 

J drr l Y L (r~~^R a )(n a - n a + Z a - Z a ) = (35) 

for all a, the potentials of these densities are zero outside their respective augmentation spheres 
(L = (l,m) is a collective angular- and magnetic quantum number). Exploiting this feature, the 
Coulomb integral reduce to 

((»» + E a Z a )) = {{n + E a za )) + Ea(O a - n a + Z a - Z a )) + 2E a (n a + Z a \n a - h a + Z a - Z a ) 
= ((n + Ea Z a )) + Ea (((»" + Z a )) - ((n a + Z a ))) 

where it has been used that inside the augmentation spheres n = n a . In this expression, we have 
circumvented all of the previous problems. None of the terms correlates functions on different grids, 
there is only a single summation over the atomic sites, and furthermore the only thing that has to 
be evaluated in the full space is the Hartree energy of n(r) + ^ a Z a (r) which is charge neutral 
(see eq. (42)). 

Inserting the final expression in (30), we see that 



(36) 



Ec[n] = i((n + E a Z a )) + \ E (((»" + za ))' + Z a )) 

a 

= u H [~p] + (((" a )) + 2 (^ a l za ) - ((" a + ^ a )) 



2 One could separate the terms in other ways, but it is impossible to separate the smooth and the localized terms 
completely. 
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where we have introduced the smooth total density 

~p(v)=n + Y,Z a {v). 



(37) 



Note that the problem with the self interaction error of the nuclei could easily be resolved once it 
was moved to the atom centered part, as handling charged densities is not a problem on radial 
grids. 

To obtain an explicit expression for the compensation charges, we make a multipole expansion 
of Z a (r) 

Z a = Y,Ql~9 a L(r), (38) 



where cj£(r) is any smooth function localized within |r — R a | < r", satisfying 



J drr l Y L (f^TL»)gl,(r) = 8 W 



(39) 



Plugging the expansion (38) into equations (35), we see that the expansion coefficients Q a h from 
must be chosen according to 



Ql = J dvr l Y L {v) K(r) - h a {r) + Z»(r)) = A a S lfi + ]T A^D? ( 



(40) 



where 



A a = dr «(r) - h a c {r)) - Z a /Vin 



a a 



drr l Y L (r)[<Pl(r)cf>l(r) - <(r)0? (r)] 



(41a) 
(41b) 



and it has been used that the core densities are spherical n"(r) = n^(r)Y 00 (i) (we consider only 
closed shell frozen cores) . This completely defines the compensation charges Z a (r) . 
Note that the special case / = of (35), implies that 



dr 



dr 



Kr) + ^Z a (r) 



= 



= / dr 



>(r) + £z°(r) 



J dr p(r) = J dr p(r) = 



(42) 



i.e. that the smooth total density has zero total charge, making the evaluation of the Hartree 
energy numerically convenient. 

In summary, we conclude that the classical coulomb interaction energy which in the KS 
formalism was given by -Ec[?i] = ^-[p], in the PAW formalism becomes a pure Hartree energy (no 
self-interaction correction) of the smooth total density p plus some one-center corrections 



E c [n] = U H {p] + J2^[{Dti 2 }] 



(43) 
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where the corrections 



&E a c [{D? li2 }] = + (n a \Z a ) - - (n a \Z a ) - \{{Z a )) 

n a c (r) 



\[{{<))-{{nt))]-Z a / dr 



1 

2 . 



k \<t>lK)-{K<l>lK)-z a / dr 



L 



V £>?* 



A a i a I i a ±a 



i 



Using that the potential of a spherical harmonic (times some radial function) is itself a spheical 
harmonic of the same angular momentum, we see that (<7£|<?£/) oc Sll 1 an d (n^\g1) oc S^q. Noting 
that Q1 by virtue of (40) is a functional of the density matrix, and inserting this, we get 



(44) 



where 



AC a =\ [(«)) - ((»«)) - (A a ) 2 ((3 a ))] - A a (« ) - V^Z a / dr 



<(r) 



A ^ 2 



,<W^ 2 (r) 



Ar" 1 — 



MO - (<«IO - Z a J dr 
A Q (CC 2 |5o a o) - Ag , lll2 (a q (« ) + ((§&))) 



E [l^W^SK) + |AJ i3i4 (<^Js2) +Al t . lfe ((s£))Ai i8i 



(45) 



(46) 



(47) 



Note that all integrals can be limited to the inside of the augmentation sphere. For example 
( ( Pi 1 ( l ) i 2 \ n c) has contributions outside the augmentation sphere, but these are exactly canceled 
by the contributions outside the spheres of \n^), in which region the two expressions are 

identical. 

The AC? „■ „■ „■ tensor has been written in a symmetric form, such that it is invariant under the 
following symmetry operations: 



ll <-» «2 



l 3 <-> l 4 



(48) 



1.5.3 Summary 

Summing up all the energy contributions, we see that the Kohn-Sham total energy 

E[n] = T s [{ip n }] + U' H [p] + E xc [n] 

can be separated into a part calculated on smooth functions, E, and some atomic corrections, AE a , 
involving quantities localized around the nuclei only. 



where the smooth part 



E = E + ^AE a 

a 

E = T s [{$ n }]+U H [p] + E xc [h] 



(49) 



(50) 
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is the usual energy functional, but evaluated on the smooth functions n and p instead of n and p, 
and with the soft compensation charges Z a instead of the nuclei charges Z a (r). The corrections 
are given by 

A£ a -AT c a + AC a + ^(AXf il2 +AC^ 2 )+ £ D% 2 AC^^D^ + AE% c [{D? li3 }] (51) 

where T", AX" i2 , ACf ii2 , and AC?^ f are system independent tensors that can be pre-calculated 
and stored for each specie in the periodic table of elements. 

Both the Hamiltonian and the forces can be derived from the total energy functional (49) as 
will be shown in the following two sections. 



1.6 The Transformed Kohn-Sham Equation 

The variational quantity in the PAW formalism is the smooth wave function ip n . From this, all other 
quantities, such as the density matrix, the soft compensation charges, the transformation operator, 
etc. are determined by various projections of onto the projector functions, and expansions 
in our chosen basis functions, the partial and smooth partial waves. To obtain the smooth wave 
functions, we need to solve the eigenvalue equation 

HMr) = e n SMr), (52) 

where the overlap operator S — T^T and H = T'HT is the transformed Hamiltonian. 



1.6.1 Orthogonality 

In the original form, the eigen states of the KS equation where orthogonal, i.e. (V'nlV'm) = 8 nm 
while in the transformed version 

(i> n \t tf $ m ) = 5 nm (53) 

i.e. the smooth wave function are only orthogonal with respect to the weight S. 
The explicit form of the overlap operator is 

S = f*f 

= (l + Ea^ , ) t (l + E Q 7' a ) 
_ i _|_ yf^ + T a + r f a ^ i j a 



= 1+ E [Eiffi>(«i - (Ci)EiC)(p? 2 i+EiC)^ 2 iE(ic> - ioki 

- (cj)E(io - \& 2 ))m] 

ii i2 

= 1 + EE WMM) - ($M))(Pt\ 

a i 1 i 2 

= 1 + EEi^)^ A So,M 2 (^ 2 i 

The orthogonality condition (53) must be kept in mind when applying numerical schemes for solving 
(52). For example plane waves are no longer orthogonal, in the sense that (Gp|G') ^ <5g.c- 



1.6.2 The Hamiltonian 

To determine the transformed Hamiltonian, one could apply the transformation H = T^HT 
directly, which would be straight forward for the local parts of H, but to take advantage of the 
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trick used to determine the total energy of the nonlocal term (Ec[n\), we make use of the relation 

6E 



f n Hip n (r). 



(55) 



Using this, we get 
SE 



5i,*(r) 6il>*(t) 

= ST s [{^ n }} 



T s [{i, n }} + E xc [h] + U H [~ P ] + AE a [{Dt ii2 }} 



dr' 



5E xc [n] | SU H [p] 
5n(r') 6n(r') 



Sh(r') 



dr 



l 5U H [n + E a Z a }SZ a (r') SAE a 



' 5Z*{v>) 5Df il2 6D? ii2 

-fn(~W Z )M*) 

+ [ dr'[v xc [n]{r') + u H [fi{r')]f n 5(r-r')$ n {r') 



EE 



SAE a 

1*2 



«(r)P, a 



where w KC [n](r) = S f*ffi is the usual local (LDA) or semilocal (GGA) exchange correlation 
potential, and ujj[n](r) = = / dr' is the usual Hartree potential. 



From these results, we can write down the transformed Hamiltonian as 
H = | V 2 + u H [p] + v xc [h] + E l^ix) M W» (P<L I > 

where the nonlocal part of the Hamiltonian is given in terms of the tensor 

SAE a 



(56) 



A K*. =E A ^ / dru H [p]{r)~gl{r) 



SDf ., 



]T &li li2 J dru H [p]{r)~gl{r) + AT^ 2 + AC^ 2 + 2 £ AC;',,, ,.ir 



SAE XC 



(57) 

Note that to justify taking the derivative with respect to D only, and not its complex conjugate, 
the symmetry properties (48) has been used to get D^AC^^D^ = D^AC^D^. 

1.7 Forces in PAW 

In the ground state, the forces on each nuclei can be calculated directly from 

dE 



pa _ 



dE 
dE 
dE 



E 



dE d\j> n ) 



+ h.c. 



(58) 



a 



E-^ e "^«i^R^i^") 



1.8 Summary 
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where h.c. denotes the hermitian conjugate. To get to the second line, the chain rule has been 
applied. The third line follows from the relation 

dE ~ ~ 

f n H\^n) = f n e n S\i> n )- (59) 



d(ip n 



The last line of (58) is obtained from the following manipulation of the orthogonality condition 
(53) 

Snm = (i>n\S\i>m) 

=> = ^(^|Sfc) = ^fc) + (^|^fc) + (^|S^ (60 ) 
* ^^ m ) + h.c. = -$ n \^ m ) 
From the expression for the overlap operator (54), it follows that 

;§ = M. (&! + <"■) («) 

which, when inserted in (58), gives the force expression 

Fa = -J^+E^«E A ^ ( P nU^n) + {i>n\^)P^) (62) 
n i t i 2 

In the case of standard xc approximations, the dependence of the total energy on the nuclei 
coordinates is 

9E = f 5E dh(r>) v dE 8Df lia f v SE d~gl(r>) 
dR a J 5n(r') dR a fr* dD? li2 dR a J ^ 8g a L {v') dR a 

giving the force expression 



(63) 



dR a 



8R a 



(64) 



1.8 Summary 

The PAW KS equation to be solved is 

H$ n ) = e n S$ n ) (65) 

with S, and H given by 

S = 1 + E E I^>^ A oVm 2 (Pfe I (66a) 
f = -iV 2 + Mff [p](r) + „ BC [»](r) +EE l^) A ^ 2 (P?J (66b) 

a i!i 2 
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where 

A ^ 2 =E A ^ 2 / ^[p](rM(r)+A^ i2 +ACf ii2 +2^AC« i2i3i4 ^ 3i4 + ^ (67) 

The total energy can then be evaluated by 

E = T s [{$ n }] + U H [p] + E xc [h] + A£ a (68) 

a 

with AE a given by 

AE a = T° + V(AZ Q , + AC? ? )Df, + V L> a * AC? , , Df , + A£±. ({£>?■ }) (69) 

21*223*4 

Having solved the eigenvalue problem (65), the eigenvalues are known. This can be used to 
determine, for example, the kinetic energy of the pseudo wave functions, T s [n], without doing the 
explicit (and computationally costly) computation. This can be seen by operating with ^ n f n {ipn\ 
on eq. (66b) to get: 

T s [{i>n}} = - / dr ^) ~ M 1 ")] [u H [p](r) +« BC [ft](r)] - J2J2 AH hiMi* ( 7 °) 



2 Implementing PAW 

For an implementation of PAW, one must specify a large number of data for each chemical element. 
This constitutes a data set which uniquely determines how the on-site PAW transformation works, 
at the site of the specific atom. For the generation of such data sets, one needs an atomic DFT 
program, with which basis sets can be generated. How to perform DFT calculations efficiently on 
an isolated atom will be discussed in the first section of this chapter, and the actual choice of data 
set parameters will be discussed in the second. The atomic DFT program plays the additional role 
of a small test program, against which implementations in the full PAW program can be tested. 



2.1 Atoms 

If we consider the Kohn-Sham equation for an isolated atom, (described by a non spin-dependent 
Hamiltonian) , it is well known that the eigenstates can be represented by the product 

^{™) = Rj{r)Y L {i) X aM) (71) 

where Rj are real radial function, and Yl are the (complex valued) spherical harmonics, i — (n, I, m), 
j = (n, I), and L — (I, m). 

Assuming identical filling of all atomic orbitals, i.e. fi a — fj, the density becomes 

i at j 

where the first factor of 2 comes from the sum over spin, and the second factor from the sum over 
the magnetic quantum number using that 

Ei^i 2 = ^ (73) 



The identical filling of degenerate states is exact for closed shell systems, and corresponds to a 
spherical averaging of the density for open shell systems. 
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Determining potentials in a spherical coordinate system is usually done by exploiting the 
expansion of the Coulomb kernel 



(74) 



with r< = min(r, r') and r> = max(r, r'). Using this it is seen that for any density with a known 
angular dependence, e.g. the density i?(r)Y£,(f), the potential can be determined by 



,R(r')Y L (v') 



v[R{r)Y L (v)]{v) = Jdv 

47T 



2l + 1 Y L (r)J Q r»*>M^ (75) 

47T 

STT y '< f > 



ii 



/r V+i Z" 00 / r \ l 
dr'R{r')r'(-) + / dr'R(r')r' 



if the angular dependence is not a spherical harmonic, one can always do a multipole expansion, 
and use the above expression on the individual terms. 

In the case of a radial density n(r) = n(r), the Hartree potential becomes 

u H {r)=— dr'n(r')r' 2 + 4tt / dr'n{r')r' (76) 

A purely radial dependent density also implies that the xc-potential is a radial function. Using 
this, the entire KS equation can be reduced to a ID problem in r, while the angular part is treated 
analytically. 

2.1.1 The Radial Kohn-Sham Equation 

For a spherical KS potential, and using that Yl are eigenstates of the Laplacian, the KS equation 
can be reduced to the simpler one-dimensional second order eigenvalue problem 



Id 2 Id 1(1+1) 



v a (r) 



Rjir) = e.R^r) (77) 



2 dr 2 r dr 2r 2 
If we introduce the radial wave function uj (r) defined by 

rR j (r) = u j (r) (78) 

the KS equation can be written as 

u'!(r) + (2e 3 - 2v s (r) - ^ ( ^) Uj .( r ) = (79) 
which is easily integrated using standard techniques. See e.g. [5, chapter 6]. 
2.2 The Atomic Data Set of PAW 

The very large degree of freedom when choosing the functions defining the PAW transformation 
means that the choice varies a great deal between different implementations. In any actual 
implementation expansions are obviously finite, and many numerical considerations must be made 
when choosing these basis sets, to ensure fast and reliable convergence. This section provides an 
overview of the information needed for uniquely defining the PAW transformation, and the level of 
freedom when choosing the parameters. 
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The Partial Waves 

The basis functions, <pf, for the expansion of \ip n ) should be chosen to ensure a fast convergence 
to the KS wave function. For this reason we choose the partial waves as the eigenstates of the 
Schrodinger equation for the isolated spin-saturated atoms. Thus the index i is a combination of 
main-, angular-, and magnetic quantum number, (n,l,m). And the explicit form is 

tf(r) = cj> a nl (r)Y lm (r) 

where (r) are the solutions of the radial KS Schrodinger equation (77) , and Y\ m are the spherical 
harmonics. For convenience we choose 4>f (r) to be real, i.e. we use the real spherical harmonics 
instead of the complex valued. This choice of partial waves implies that the smooth partial waves 
and the smooth projector functions can also be chosen real, and as products of some radial functions 
and the same real spherical harmonic. 

Note that including unbound states of the radial KS equation in the partial waves is not a 
problem, as the diverging tail is exactly canceled by the smooth partial waves. In practice we only 
integrate the KS equation outward from the origin to the cutoff radius for unbound states, thus 
making the energies free parameters. In principle the same could be done for the bound states, but 
in GPAW, the energies of bound states are fixed by making the inward integration for these states 
and doing the usual matching (see e.g. [5, chapter 6]), i.e. the energies are chosen as the eigen 
energies of the system. 



The Smooth Partial Waves 

The smooth partial waves ipf (r) are per construction identical to the partial waves outside the 
augmentation sphere. Inside the spheres, we can choose them as any smooth continuation. Presently 
GPAW uses simple 6'th order polynomials of even powers only (as odd powers in r results in a kink 
in the functions at the origin, i.e. that the first derivatives are not defined at this point), where 
the coefficients are used to match the partial waves smoothly at r — r c . Other codes uses Bessel 
functions [4] or Gaussians. 



The Smooth Projector Functions 

The smooth projector functions are a bit more tricky. Making them orthonormal to <fif(r) is a 
simple task of applying an orthonormalization procedure. This is the only formal requirement, but 
in any actual implementation all expansions are necessarily finite, and we therefore want them to 
converge as fast as possible, so only a few terms needs to be evaluated. 

A popular choice is to determine the smooth projector functions according to 



(80) 



or equivalcntly 



2dr 2 



1 d 

r dr 



1(1+1) 
2r 2 



v s (r) 



(81) 



where v s (r) is the smooth KS potential ujj[p](r) + v xc [h](r). And then enforce the complementary 



orthogonality condition (p^W,) 



inside the augmentation sphere, e.g. by a standard Gram- 



Schmidt procedure. Using this procedure ensures that the reference atom is described correctly 
despite the finite number of projectors. 



The Smooth Compensation Charge Expansion Functions 

The smooth compensation charges g£(r), are products of spherical harmonics, and radial functions 
gf{r) satisfying that 

r drr l Y L (r)g a L ,(r) = S LL , (82) 
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In GPAW the radial functions are chosen as generalized Gaussian according to (here shown for 
R Q = 0): 

~gl(r) = ~gt(r)Y L (r) , §f(r) = ^^^(4« Q )'+ 3 /V e -^ 2 (83) 

where the atom-dependent decay factor a is chosen such that the charges are localized within the 
augmentation sphere. 

The Core- and Smooth Core Densities 

The core density follows directly from the all electron partial waves by 

core core 

n c (r) = ]T |<Mr)| 2 = £2(2^ + l)|^(r)| 2 /4^ (84) 

i 3 

The smooth core densities n"(r) are like the smooth partial waves expanded in a few (two or 
three) Bessel functions, Gaussians, polynomials or otherwise, fitted such that it matches the true 
core density smoothly at the cut-off radius. 

The Localized Potential 

An additional freedom in PAW is that for any operator L, localized within the augmentation 
spheres, we can exploit the identity (8) 

£l#)(p?l = i ( 85 ) 

i 

valid within the spheres, to get 

a i 1 i 2 

so for any potential v(r) = ^ a v a (r — R a ) localized within the augmentation spheres (i.e. v a (r) = 
for r > r") we get the identity 



= y dm(r)^tj a (r)-^ J drh a v a (r) 



This expression can be used as an 'intelligent zero'. Using this, we can make the replacement of 
the smooth potential 

v s (r) = u H [p](r) + v xc [h)(v) v s (v) = u H [p](r) + v xc [n)(r) + v(r) (86) 

if we at the same time add 

B a + V ABf , Df , (87) 

to the energy corrections AE a , where 

B a = -J dvn a c v a {v) and AB^ 2 =-J dr<^« a (r) 

This also implies that Bfi should be added to AHf ii2 . 

The advantage of doing this is that the Hartree potential and the xc-potential might not be 
optimally smooth close to the nuclei, but if we define the localized potential properly, the sum of 
the three potentials might still be smooth. Thus one can initially evaluate uh[p\{y) and Ua;c[^]( r ) 
on an extra fine grid, add v(r) and then restrict the total potential to the coarse grid again before 
solving the KS equation. 

The typical way of constructing the localized potentials v a is by expanding it in some basis, 
and then choosing the coefficients such that the potential ti^[p](r) + w :EC [n](r) + v(r) is optimally 
smooth at the core for the reference system. 

Inclusion of v a (r) changes the forces on each atom only through the redefinitions of v s (r) and 
AH a ■ . 

UI2 
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Summary 

When constructing a data set for a specific atom, one must specify the following quantities, all 
defined within the augmentation spheres only: 

1. (f>f from radial KS equation 

2. (j)® by appropriate smooth continuation of <pf 

3. p$ from equation (80) 

4. gl localized within r < r c , and satisfying J drr 1 (/£(r)Yf/(r — R a ) = 6ll> 

5. n" follows from <j>? by (84) 

6. n a c by appropriate smooth continuation of fi" 

7. v a localized within r < r", otherwise freely chosen to make v s optimally smooth for the 
reference system 

The adjustable parameters besides the shape of 4> a , 9li ano - are 

1. Cut-off radii r" (which can also depend on i) 

2. Frozen core states 

3. Number of basis set functions (range of index i on 0°, <^°, and pf) 

4. Energies of unbound partial waves 

Choosing these parameters in such a way that the basis is sufficient for the description of all 
possible environments for the specific chemical element, while still ensuring a smooth pseudo wave 
function is a delicate procedure, although the optimal parameter choice is more stable than for e.g. 
norm conserving or ultra soft pseudopotentials. 

Once the quantities above have been constructed, all other ingredients of the PAW transformation 
follows, such as A\ A^„ T«, AT? lh , AC°, AC» !2 , AC? Mi , Av a , and Av£ ia . The first two 
are needed for the construction of the compensation charges and the overlap operator, and the rest 
for determining the Hamiltonian, and for evaluating the total energy. 

3 Non-standard Quantities 

The preceding sections have described the details of making a standard DFT scheme work within 
the PAW formalism. This section will focus on what the PAW transform does to quantities needed 
for post-processing or expansions to DFT. 

It is a big advantage of the PAW method, that it is formally exactly equivalent to all-electron 
methods (with a frozen core) but is computationally comparable to doing pseudopotential cal- 
culations. In pseudopotential approaches, projecting out the core region is handled by a static 
projection kernel, while in PAW this projection kernel is dynamically updated during the SCF-cycle 
via an expansion of the core region in a local atomic basis set. This has the drawback for the end 
user, the equations for all quantities most be modified to account the dual basis set description. 

3.1 The External Potential in PAW 

As an example of the principle in accommodating expressions to the PAW formalism, we will here 
consider the application of an external potential in DFT. 

Without the PAW transformation, this addition is trivial, as the desired potential, u cxt (r), 
should simply be added to the effective KS potential, and the total energy adjusted by the energy 
associated with the external potential E ext = J drn(r)u cxt (r). 



3.2 All-electron Density 



19 



In PAW, the density decomposes into pseudo and atomic parts, so that 

#ext= / drn(r)v cxt (r) +J2 / * [n a (r) - n a (r)] v cxt (r). 

Implying that both a pseudo energy contribution E cxt — J drh(r)v ex t(r) and atomic corrections 
AE£ xt = J dr [n a (r) — h a (r)] ?w(r) should be added to the total energy. 
In PAW, the Hamiltonian has the structure: 



H 



1 



dE 



fn\lpn) d(lp n \ 



a i\i2 

In our case, the extra contributions due to the external potential are: 

ffext(r) = fcxt(r) 

and 

AH^f = J dr« ext (r){Ci(r)^ 2 (r)-C(r)^ 2 (r)} 
Thus we can write the atomic energy contribution as: 



(89) 



drv cxt (r) 



i«(r) - n«(r) + ]>>U - ^(r)0? 2 (r)} 



= J drv cxt (r) K(r) - n a c (r)} + ^ D^AH^ 



To evaluate the first term in the last line, the external potential should be expanded in some radial 
function at each nuclei e.g. the gaussians g£, as the integral of these with the core densities is 
already precalculated. 

For example, a zero-order (monopole) expansion, equivalent to the assumption 



w«t(r) «f cx t(R a ) , for |r-R a | < r a c 



Leads to the expression: 



AE« xt = v cxt (R a )(V^Q a ao + Z a ) 



ah c : 



Ucxt (R a )V47rAg 0jn 



Linear external potentials (corresponding to a homogeneous applied electric field) can be handled 
exactly by doing an expansion to first order. This has been used in GPAW in e.g. the paper [6] . 

3.2 All-electron Density 

During the self-consistency cycle of DFT, the all-electron quantities are at all times available 
in principle. In practise, they are never handled directly, but rather in the composite basis 
representation of a global pseudo description augmented by local atomic basis functions. For some 
post processing purposes it can however be desirable to reconstruct all-electron quantities on a 
single regular grid. 

As an example, consider the all-electron density, which can formally be reconstructed by 



i(r) = n(r)+Y^ 



^(r) + ]T D? ii2 {4>l (r)^ 2 (r) - ^ (r)<^ (r) 



Transferring this to a uniform grid will of coarse re-introduce the problem of describing sharply 
peaked atomic orbitals on a uniform grid, but as it is only needed for post processing, and not in 
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the self-consistency, it can be afforded to interpolating the pseudo density to an extra fine grid, 
before adding the summed atomic corrections from the radial grid. 

One common use of the all-electron density is for the application of Bader analysis [7]. The 
advantage of applying this to the all-electron density instead of the pseudo density, is that it can 
be proved that the total electron density only has maxima's at the nuclei, such that there will only 
be one Bader volume associated with each atom. This does not hold for the pseudo density, which 
can result in detached Bader volumes. In addition, the dividing surfaces found if applied to the 
pseudo density will be wrong if these intersect the augmentation sphere. 

In practice, the reconstructed total density will not integrate correctly due to the inaccurate 
description of a uniform grid in the core regions of especially heavy elements. But as the numerically 
exact value of the integral over the atomic corrections are known from the radial grid description 
(= 4:irJ2ij D?jA1 A j ), this deficiency can easily be remedied. As long as the density is correctly 
described at the dividing surfaces, these will still be determined correctly. 

3.3 Wannier Orbitals 

When constructing Wannier functions, the only quantities that needs to be supplied by the DFT 
calculator are the integrals Z^ = (ip ni |e~ G ' r |-0n 2 )j where G is one of at most 6 possible (3 in an 
orthorhombic cell) vectors connecting nearest neighbor cells in the reciprocal lattice. [8, 9] 
When introducing the PAW transformation, this quantity can be expressed as 

ZZ ina = $m\e- G - T \4>n 2 ) +^E F "nC 2 (<C|e- G ' r |C> " (< l^'IC)) ■ 

a iii 2 

Even for small systems, the phase of the exponential of the last integral does not vary significantly 
over the region of space, where pf is non-zero. The integral in the last term can therefore safely be 
approximated by 

-GR a pa* pa /TlAa 

iii 2 

3.4 Local Properties 

This section describes quantities that can somehow be related to a specific atom. As the PAW 
transform utilizes an inherent partitioning of space into atomic regions, such quantities are usually 
readily extractable from already determined atomic attributes, such as the density matrices or the 
projector overlaps P^, which are by construction simultaneous expansion coefficients of both the 
pseudo and the all-electron wave functions inside the augmentation spheres. 

3.4.1 Projected Density of States 

The projection of the all electron wave functions onto the all electron partial waves, (i.e. the all 
electron wave functions of the isolated atoms) <fif , is within the PAW formalism given by 

(<f>Wn) = (^li)+^E^I^((<l<> - (k'M))(Pt\^n) (90) 
a' i±i 2 

Using that projectors and pseudo partial waves form a complete basis within the augmentation 
spheres, this can be re-expressed as 

(m n ) = + £ ^2(&w[)&s£ i2 p< 2 (9i) 

a'^a i\i 2 

if the chosen orbital index 'i' correspond to a bound state, the overlaps ((j>i\pf), a! ^ a will be 
small, and we see that we can approximate 



(92) 



3.5 Coulomb Integrals 
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The coefficients = (pt\ip n ), can thus be used as a qualitative measure of the local character of 
the true all electron wave functions. As the coefficients are already calculated and used in the SCF 
cycle, it involves no extra computational cost to determine quantities related directly to these. 
These can be used to define an atomic orbital projected density of states 

n 4 ( £ ) = ^<5( £ -6„)|P„ a J 2 . (93) 



3.4.2 Local Magnetic Moments 

As the projection coefficients are simultaneous expansion coefficients of the pseudo and the all- 
electron wave functions inside the augmentation spheres, it can be seen that inside these, the 
all-electron density is given by (for a complete set of partial waves) 

n(r)=^Df li3 <(r)C 2 (r)+<(r), |r-R a |<r c Q . (94) 
This can be used to assign a local magnetic moment to each atom according to 

where AN is an integration over products of AE waves truncated to the interior of the augmentation 
sphere 

AAf il2 = / dr^(r)^(r). 

Note that this will not add up to the total magnetic moment J dr(rif(r) — nj^(r)), due to the 
interstitial space between augmentation spheres, and must be scaled if this is desired. 



3.4.3 LDA + U 

The atom projected density matrix Df ^ can also be used to do LDA -I- U calculations. The GPAW 
implementation follows the LDA + U implementation in VASP[10], which is based on the particular 
branch of LDA + U suggested by Dudarev et aZ. [11], where you set the effective (U-J) parameter. 
The key notion is that from (94) one can define an (valence-) orbital density matrix 

$U = IC)A ? ii 2 vA ? J- 

Thus doing LDA + U is a simple matter of picking out the d-type elements of D a , and adding to 
the total energy the contribution 

d typo / \ 

E E 2^ { D ^-Y, D lisKiA (95) 

and adding the gradient of this to the Hamiltonian 

d typo j j 

E E i#i>2 ( s ^- 2D Zi*)<pi\ ( 96 ) 

a iii2 



3.5 Coulomb Integrals 

When trying to describe electron interactions beyond the level of standard (semi-) local density 
approximations, one will often need Coulomb matrix elements of the type 
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where the orbital pair density n nn /(r) — ?/>*(r)^v(r). 

Such elements are needed in some formulations of vdW functionals (although not the one 
implemented in GPAw), in linear-response TDDFT (see e.g. [12]) where only pair densities 
corresponding to electron-hole pairs are needed, in exact exchange or hybrid functionals (see next 
section) where only elements of the form K nn i )nn i where both indices correspond to occupied states, 
are needed, and for GW calculations (see e.g. [13]), where all elements are needed. 

Introducing the PAW transformation in (97), the pair densities partition according to 

rw(r) = rw(r) + £ « n ,(r) - f»° n ,(r)) (98) 

a 

with the obvious definitions 

h nn ,=rj n , c=Ec;^ 2 cc c-=Ec;^ 2 cc- (99) 

Exactly like with the Hartree potential, direct insertion of this in (97) would, due to the non-local 
nature of the Coulomb kernel, lead to undesired cross terms between different augmentation spheres. 
As before, such terms can be avoided by introducing some compensation charges, Z* n t, chosen such 
that the potential of nf ln , — ii° n , — Z® n , are zero outside their respective augmentation spheres. 
This is achieved by doing a multipole expansion and requiring the expansion coefficients to be zero, 
and entails a compensation of the form 

K,Ar) = Ql,nn~9l{r), Ql, nn , = J2 ^LMh P nt P nH 2 (100) 

(the constants A 2 ili2 are identical to those in (41b)). 

Introduction of such compensation charges makes it possible to obtain the clean partitioning 

K nn ' <mm i — (Pnn'lPmm 1 ) + 2 ^""^ P mii P ni 2 ^^iii 2 i3U P n'i 3 P m'ii 1 (101) 

a 1112*3*4 

Here the last term is a trivial functional of the expansion coefficients involving only the 
constants AC? i2i3i4 already precalculated for the atomic corrections to the Coulomb energy (47). 
The only computationally demanding term relates to the Coulomb matrix clement of the smooth 
compensated pair densities pij — hij +J2 a ^ij> which are expressible on coarse grids. 

The formally exact partitioning (101) makes it possible, at moderate computational effort, to 
obtain Coulomb matrix elements in a representation approaching the infinite basis set limit. In 
standard implementations, such elements are usually only available in atomic basis sets, where the 
convergence of the basis is problematic. At the same time, all information on the nodal structure of 
the all-electron wave functions in the core region is retained, which is important due the non-local 
probing of the Coulomb operator. In standard pseudopotential schemes, this information is lost, 
leading to an uncontrolled approximation to K nn i mm i. 

As a technical issue, we note that integration over the the Coulomb kernel l/|r — r'| is done by 
solving the associated Poisson equation, as for the Hartree potential, whereby the calculation of each 
element can be efficiently parallelized using domain decomposition. The integral J drp nn > (r) = 5 nn > 
shows that the compensated pair densities p nn have a non-zero total charge, which is problematic 
for the determination of the associated potential. For periodic systems, charge neutrality is enforced 
by subtracting a homogeneous background charge, and the error so introduced is removed to 
leading order (V^ 1 ^ 3 where V the the volume of the simulation box) by adding the potential of a 
missing probe charge in an otherwise periodically repeated array of probe charges embedded in a 
compensating homogeneous background charge. This can be determined using the standard Ewald 
technique, and corresponds to a rigid shift of the potential. For non-periodic systems, all charge is 
localized in the box, and the Poisson equation can be solved by adjusting the boundary values 
according to a multipole expansion of the pair density with respect to the center of the simulation 
box. A monopole correction is correct to the same order as the correction for periodic cells, but 
the prefactor on the error is much smaller, and leads to converged potentials even for small cells. 
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3.5.1 Exact Exchange 

The EXX energy functional is given by 



~2 ^ ^ fnfm&<T n ,<T m -K nnl nnl . (102) 



inn 



Terms where n and m both refer to valence states transform in PAW as in equation (101). Terms 
where either index refers to a core orbital can be reduce to trivial functionals of P^, resulting in 
(see e.g. [14]) 

^ val 



nin 



E 



(103) 



The term involving the AC a tensor is the PAW correction for the valence- valence interaction, and 
is similar to the correction in the equivalent expression for the Hartree energy, except that the 
order of the indices on the density matrices are interchanged. The term involving the X a tensor 
represents the valence-core exchange interaction energy. E££~ c is simply the (constant) exchange 
energy of the core electrons. 

The system independent Hermitian tensor Xf iit> is given by: 

*«, J J ^ jT^j 

core" / X ff , (104) 

= EE aifl (E G ^ G ^) J J drdr'^ul(r)ul(r)ul(r')ul(r'). 

Although the valence-core interaction is computationally trivial to include, it is not unimportant, 
giving rise to shifts in the valence eigenvalues of up to leV (though only a few kcal/mol in 
atomization energies), and we note that this contribution is unavailable in pseudopotential schemes. 
The core-core exchange is simply a reference energy, and will not affect self-consistency or energy 
differences. 

For the iterative minimization schemes used in real-space and plane wave codes, the explicit 
form of the non-local Fock operator w NL (r,r') is never needed, and would indeed be impossible 
to represent on any realistic grid. Instead only the action of the operator on a state is needed. 
As with the Hamiltonian operator, the action on the pseudo waves is derived via the relation 
/ n 6 NL |VVi) = dE^/d(tp n \. Referring to [14] for a derivation, we merely state the result 

u NL |V>n) = E f™v nm (*)\^ m ) 



EEift) 

a iii 2 
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(105) 



is the solution of V 2 w„ m (r) = -4vrp nm (r), and v% mMia = J2l A h ll2 1 dvg a L {v)i nm {v) . 
Again the computationally demanding first term is related to smooth pseudo quantities only, 
which can be accurately represented on coarse grids, making it possible to do basis set converged 
self-consistent EXX calculations at a relatively modest cost. Applying the Fock operator is however 
still expensive, as a Poisson equation must be solved for all pairs of orbit als. 

As a technical consideration, note that the effect of the atomic corrections due to valence- valence, 
valence-core, and core-core exchange interactions can simply be incorporated into the standard 
equations by redefining equations (47), (46), and (45) respectively, which will also take care of 
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the last two terms in the Fock operator above. The introduction of the pair orbital compensation 
charges does however lead to a non-trivial correction to the Fock operato; the term proportional to 
v nm % % • This term also leads to a distinct contribution when calculating the kinetic energy via 
the eigenvalues as done in equation (70). The additional term (besides those related to redefining 
(45)-(47)) 



/m<Wm / dr{w(r)^(r)^ m (r) - EE KnPmi 2 <m,ni 



(106) 



should be added to the right hand side of (70) on inclusion of exact exchange. 

In a similar fashion, the compensation charges leads to an additional force contribution in 
equation (64) given by 



^xx — y ] fnfnm$a n a m ) / ^ r ^™m( r ) ^ ] ^ > nh^'mi2 ^ ] ^Li t i 



d~gt(v') 



(107) 
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3.5.2 Optimized Effective Potential 

The optimized effective potential (OEP) method, is a way of converting the non-local Fock operator 
into a local form — u£(r). 

One way to derive the OEP equations in standard KS-DFT, is to use perturbation theory along 
the adiabatic connection (Gorling-Levy perturbation theory [15]). 

On converting the OEP equation to the PAW formalism, it should be remembered that local 
potentials in PAW transform to a local pseudo part plus non-local atomic corrections. Hence we 
want to arrive at a potential of the form 

vl = ^(r)+E5>£> A <^"J< ( 108 ) 

a i 1 i 2 

where both the pseudo part u£ as well as the coefficients Auf i should be determined. 

The derivation is more or less straight forward, if one remembers the the PAW KS equation is 
a generalized eigenvalue problem, that the variational quantity is the pseudo orbitals, and that the 
first order shift in the density has both a pseudo and an atomic part. The result is 

£ /„C(r) £ ^(r ^^J^ + c.c. = (109a) 
E « E <* m| f 7*^> + c.c. = (109b) 

n m^tn 

where is the non-local exchange operator of equation (105) and is the local version in (108). 

These can be solved iteratively starting from a local density-function approximation to the 
exchange potential in the spirit of [16]. 

It might seem that OEP is just extra work on top of the already expensive non-local operator, 
but it can in some cases be faster, as the number of SCF iterations in the KS cycle are greatly 
reduced. 
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